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Executive  Summary: 


Research  is  conducted  to  find  a  robust  wavelet 
based  algorithm  for  automatic  identification  of 
push-to-talk  radio  transmitters.  Digital  data  at  an 
intermediate  frequency  (IF)  is  pre-processed  to 
translate  it  into  a  form  applicable  to  wavelet 
analysis.  The  magnitudes  and  the  relative  positions  of 
the  extrema  from  the  best  suited  wavelet  scales  are 
used  in  conjunction  with  a  distance  measure  algorithm 
to  determine  from  which  transmitter  any  particular 
signal  may  have  originated.  The  distance  measure 
algorithm  can  correctly  identify  the  four  signal  sets 
provided  for  the  research.  Even  at  reduced  signal  to 
noise  (SNR)  ratios  good  identification  is  obtained. 
Automatic  selection  of  a  good  (i.e.  optimal)  threshold, 
the  simultaneous  use  of  several  scale  outputs  and  a 
more  robust  choice  of  the  reference  templates  should  be 
addressed  in  follow  on  work.  The  current  version  of  the 
algorithm  can  classify  all  sets,  with  at  least  an  80 
percent  reliability  at  an  SNR  on  the  average  10  dB 
worse  than  the  original  recordings.  The  procedure  lends 
itself  to  be  automated,  minimizing  human  interaction. 
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I .  Introduction 

Signals  generated  by  different  transmitters  are  believed  to 
possess  somewhat  different  characteristics.  The  transient 
response  of  a  transmitter  is  defined  to  be  the  carrier  strength 
change  from  the  off- state  to  the  on- state,  and  of  course  from  the 
on- state  to  the  off -state.  The  turn  on  response  is  unique  to  a 
specific  transmitter  and  may  even  differ  among  units  of  the  same 
make  and  model.  These  observations  are  not  applicable  to  the 
turn-off  response.  Therefore,  the  uniqueness  of  the  tum-on 
transients  is  exploited  to  identify  the  source  of  the  signals  and 
therefore  the  transmitter. 

Time -frequency  analysis  of  stationary  signals  is  a  well 
studied  and  known  subject.  The  Fourier  Transform  (FT)  method  is 
well  suited  for  this  type  of  analysis.  However,  when  the  signal 
of  interest  is  non- stationary,  the  FT  method  is  not  appropriate, 
since  it  uses  a  complex  exponential  basis  function  that  exists 
over  infinite  time.  A  sliding  time-window  (Gaussian)  was 
introduced  by  Gabor  (1946)  to  gain  time  information  from  the  FT 
method.  This  modified  Fourier  Transform  method  is  called  the 
Short  Time  Fourier  Transform  (STFT)  .  Once  a  window  is  chosen,  the 
time- frequency  resolution  is  fixed.  This  method  requires  that  the 
signal  is  stationary  within  these  intervals. 

A  technique  more  suitable  for  transient  signal  analysis  is 
the  Wavelet  Transform  (WT)  method.  It  is  more  revealing  than  the 
STFT  method  in  terms  of  time  and  frequency  information.  Basis 
functions  of  the  WT,  unlike  the  complex  exponential  of  the  FT,  or 
STFT  are  shorter  in  time  duration  then  the  analysis  interval. 

This  compact  support  makes  the  WT  localized,  in  frequency  and  in 
time.  Moreover,  wavelets  provide  the  flexibility  to  choose  the 
particular  wavelet  function  that  is  appropriate  for  a  specific 
application.  This  is  possible  since  there  are  a  large  number  of 
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compactly  supported  wavelets  that  can  be  used  as  orthogonal  basis 
functions . 

The  purpose  of  this  report  is  to  investigate  the  use  of  the 
WT  method  to  classify  transmitter  signatures. 


II.  Signal  Preprocessing  and  Filtering. 

1.  Signal  Preprocessing: 

Representative  data  from  each  source,  that  is  one  signal 
from  each  transmitter  is  shown  in  Figure  1  and  2.  Figure  1 
displays  typical  turn  on  transients,  while  Figure  2  shows  the 
respective  turn  off  transients.  The  data  is  recorded  after  it  is 
intercepted  by  the  antenna  and  processed  by  a  typical  radio 
receiver.  The  carrier  frequency  of  the  radio  is  138.525  MHz.  The 
signals  are  filtered  with  a  1  Mhz  bandwidth  (BW) ,  digitized  at  a 
sampling  frequency  of  5  Mhz  at  a  center  frequency  of  1.075  Mhz. 
The  signals  are  in  binary  form  with  10  bits  available  for  a 
discrete  value  representation.  These  signals  (i.e.  digital 
recordings)  are  pre-processed  (i.e.  modified)  to  change  the  data 
into  a  form  suitable  for  WT  processing. 

There  are  four  steps  in  the  pre-processing  phase, 
sequentially  given  by: 

1)  taking  the  envelope, 

2)  filtering, 

3)  differentiating,  and 

4)  a  final  filtering. 

Prior  to  taking  the  envelopes,  the  DC  terms  are  removed. 

100 -point  boxcar  averagers  are  used  for  filtering  the  envelope, 
whereas  50 -point  boxcar  averagers  are  used  after  the 
differentiation.  The  sizes  of  the  filters  were  experimentally 
determined. 
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The  processing  of  the  turn  off  transient  did  not  lead  to 
identification,  hence  our  work  deals  only  with  the  turn  on 
transient  portions  of  the  data.  Widening  the  processing  band 
width  by  using  higher  sanpling  rates  may  allow  classification  of 
the  transmitters  based  on  turn  off  characteristics.  Of  course,  a 
larger  bandwidth  allows  more  noise  to  pass,  making  this  approach 
less  reliable  at  a  given  SNR  level.  The  bandwidth  (about  2.5  MHz) 
of  the  data  set  may  not  permit  the  observation  of  very  short 
duration  transients  if  they  indeed  do  exist. 


0  2000  4000  6000  0  2000  4000  6000 

k  k 


(C)  (d) 


Figure  1.  Turn  On  Transients  From  Four  Transmitters. 
Transmitters:  (a)  =  1, (b)  =2,  (c)  =3,  (d)  =4  . 
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Figure  2 .  Turn  Off  Transients 
Transmitters:  (a)  =  1, (b) 


From  Four  Transmitters. 
=  2,  (c)  =3,  (d)  =  4  . 


A  closer  look  at  the  turn  on  transients  reveals  a  signal 
with  a  typical  ramp  type  behavior  (i.e.  unit  step  response 
of  a  second  or  higher  order  linear  system)  .  They  may  differ  in 
slope,  have  unique  dips  or  a  slow  oscillation  of  the  envelope 
(i.e.  1  to  2  cycles  over  the  ramp  duration)  .  However  none  of 
these  transients,  even  with  the  DC  component  removed,  is  in  a 
form  to  which  the  WT  is  directly  applicable.  Wavelet  Transforms 
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are  not  useful  in  analyzing  low  frequency  signals,  but  they  are 
well -suited  for  short  duration  phenomena.  Thus,  it  was  necessary 
to  transform  the  data  into  a  form  suitable  for  wavelet  analysis. 
Details  are  presented  in  the  next  chapter.  It  should  be  noted 
that  a  denoising  process  tends  to  enhance  the  identification 
performance . 

To  transform  the  envelopes  of  the  signals  into  a  form  which 
allows  successful  WT  analysis,  the  envelope  is  differentiated 
changing  the  ramps  (i.e.  step  responses)  into  pulse  like  signals. 
Figure  3  shows  the  typical  results  after  the  pre-processing 
operation  when  applied  to  the  signals  shown  in  Figure  1.  The 
differentiation  is  a  high  pass  operation  that  provides  a  gain 
that  is  linearly  proportional  to  frequency.  That  means,  the  high 
frequency  components  become  emphasized.  Unfortunately,  this  also 
applies  to  the  additive  noise  components  so  that  careful 
filtering  must  be  employed  during  the  stages  of  the  pre¬ 
processing.  From  empirical  studies  and  some  theoretical 
considerations  a  filtering  operation  is  selected  and  applied  to 
the  data  before  and  after  the  differentiation. 

The  data  used  in  the  analysis  was  collected  and  recorded  by 
the  Naval  Security  Group  Activity,  Charleston,  SC.  Nine 
tum-on/off  samples  of  each  of  four  transmitters  were  recorded. 
All  the  radios  are  Motorola  models.  Each  radio  is  identified  by 
its  model  name  and  number  and  is  listed  below: 


Radio: 

Transmitter 

1 

(Trl) 

Model : 

Maxtrac 

Transmitter 

2 

(Tr2) 

Saber 

Transmitter 

3 

(Tr3) 

HT440 

Transmitter 

4 

(Tr4) 

Saber. 

Tr2  and  Tr4  are  different  radios  of  the  same  model. 
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Figure  3.  Pre-processed  Turn  On  Signals  (  4  Transmitters) 
Transmitters:  (a)  =  1, (b)  =  2,  (c)  =  3,  (d)  =  4  . 


2 .  FILTERS 

In  the  implementation  of  the  pre-processing  two  filtering 
operations  are  performed.  The  first  one  follows  the  envelope 


operation  and  removes  some  of  the  broadband  noise.  This  filter 
(boxcar  averager  of  length  100)  is  followed  by  a  differentiation 
operation  (i.e.  first  order  difference) .  The  differentiation 
emphasizes  the  high  frequency  noise  at  the  same  rate  as  the  first 
filter  tends  to  attenuate  it.  Hence  a  second  filter  (low  pass) 
operation  is  implemented.  The  second  filter  uses  again  a  box  car 
averager,  in  this  case  of  length  50.  Boxcar  averagers  are  also 
called  integrate  -  and  -  dump  filters.  Sometimes  they  are 
referred  to  as  moving  average  filters.  They  slide  along  the  time 
series  (the  data)  and  form  an  average  over  the  length  of  data 
span  contained  within  the  filter  memory.  The  frequency  response 
of  the  boxcar  averager  is  a  digital  sine  function,  which  is  very 
similar  to  the  familiar  analog  sine  function.  One  could 
incorporate  all  three  operations  (lowpass  -  differentiation  - 
lowpass)  into  one  sophisticated  filter,  but  for  ease  of 
computation  and  efficiency  this  is  not  done.  Processing  that  uses 
median  filtering  in  place  of  the  last  lowpass  filter  in  the  pre¬ 
processing  phase  is  an  alternative  scheme  that  can  be  pursued  at 
some  time.  This  would  minimize  the  spreading  due  to  the 
convolutional  filtering  (median  filtering  tends  to  preserve  the 
rise  and  fall  times  of  a  pulse)  . 


III.  WAVELET  TRANSFORM  (WT) 

1 .  Introduction 

Wavelets  are  proportional  bandpass  processing  schemes  which 
are  also  known  as  constant  Q-filtering  (Akansu  and  Haddad,  1992; 
Burrus  and  Gopinath,  1993)  .  Typical  basis  functions  are  the 
Walsh,  Daubechies,  spline,  and  sine  basis  function.  The  Walsh 
function  is  a  set  of  rectangular  basis  function,  the  Daubechies 
functions  are  solutions  to  the  scaling  function,  the  spline 
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function  uses  a  triangular  scaling  coefficients,  comparable  to  a 
Bartlett  weighting  in  a  FIR  low  pass  filter.  The  sine  function  is 
equivalent  to  a  brick  (ideal)  low  pass  filter.  The  bandpass 
filter  structure  uses  these  basis  functions  as  impulse  response 
in  proto-type  FIR  filters  to  achieve  the  proportional  band  pass 
filtering. 

The  Wavelet  Transform  (WT)  is  founded  on  a  set  of  specific 
basis  functions,  which  are  called  wavelets  (Young,  1993)  .  They 
include  short  duration/high  frequency  and  long  duration/low 
frequency  functions.  Each  element  in  the  wavelet  set  is 
constructed  from  the  same  function,  which  is  called  the 
'analyzing  wavelet'  or  the  'mother  wavelet'. 

There  are  three  conditions  for  a  function  to  be  a  mother 
wavelet.  It  must  be  oscillatory,  it  must  decay  to  zero,  and  it 
must  integrate  to  zero  (Young,  1993) . 

The  processing  scheme  is  adopted  from  Newland,  which 
requires  that  the  scaling  coefficients  have  to  be  given  to  the 
wavelet  transform  routine  (Pitta,  1995) .  Since  the  signals  of 
interest  are  expected  to  have  a  some  discontinuity  in  phase 
and/or  in  amplitude,  basis  functions  are  selected  with  matched 
filtering  as  motivation. 

Several  WT  were  tried  with  the  Daubechies  polynomial  of 
order  8  best  suited  for  the  data  at  hand  (Payal,  1995) .  This  is 
understandable  in  a  matched  filter  sense,  if  we  compare  a  typical 
output  waveform  from  the  pre-processing  procedure  with  the  family 
of  Daubechies  polynomials.  The  processing  is  accomplished  by 
forming  the  inner  product  of  the  scaling  coefficients  (i.e.  low 
pass  filter)  and  the  wavelet  coefficients  (i.e.  high  pass  filter) 
with  the  data  as  indicated  by  equation  3.1: 
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dmn  =<  X(t) 

tM(t)  =  ao~m/2  *  (ao“e  -  nb0) 

<4o*«  (« 

(3.1) 

i|r(t)  =  2  J^d(n)  <|) (2 t~n) 

<|>(t)  =  2  J^c(.n)  <|) ( 2 t—n) 
d{n)  =  (-l)n  c(iV-n); 

where  c(n)  &  d(n)  are  the  weights  selected  for  the  lowpass  & 
highpass  filter,  respectively.  N  is  the  length  of  the  filters  and 
d^n  is  the  wavelet  transform  coefficient  at  scale  m  and  time 
(delay)  n.  The  d(n)  coefficients  are  the  wavelet  coefficients 
which  related  to  the  scaling  coefficients  in  a  sirrple  manner 
(i.e.  position  reversed  and  alternating  in  sign)  .  T  (t)  and  $  (t) 
are  the  wavelet  function  and  scaling  function,  respectively. 


2.  Filter  Banks  and  Discrete  Wavelet  Transform 

Multiresolution  analysis  can  be  implemented  by  using  a 
technique  called  Multiresolution  Pyramid  Decomposition  or 
Mallat's  algorithm  (Vetterli  and  Kovacevic,  1995).  Mallat's 
Pyramid  Algorithm  is  used  to  obtain  the  Discrete  Wavelet 
Transform  (DWT)  .  The  Discrete  Wavelet  Transform  coefficients  at 
scale  j  are  obtained  by  convolving  the  coefficients  at  scale  j+1 
with  ho(n)  and  hx(n)  .  The  impulse  responses  ho(n)  and  hx  (n)  are 
the  time  reversed  coefficients  of  c(n)  and  d(n) ,  respectively. 
This  followed  by  a  decimation  procedure  to  produce  the  expansion 
coefficients  at  scale  j .  Figure  4  shows  the  implementation  of  the 
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Mallat  Pyramid  Algorithm  for  three  levels  (i.e.  scales)  .  The 
notation  LP  represents  a  lowpass  FIR  filter,  while  HP  represents 
a  highpass  FIR  filter. 


Figure  4.  Three  Levels  of  Multiresolution  Analysis 

Due  to  its  pipe  line  structure  the  algorithm  is  very  fast. 
It  can  be  faster  than  an  FFT,  since  the  processing  cost  is 
linearly  proportional  to  the  data  length.  All  filters  in  a  given 
class  (i.e.  LP  or  HP)  are  identical,  making  the  implementation 
fairly  simple. 


IV.  IDENTIFICATION  PHASE 

1.  Reduced  Set  Representation 

One  of  the  main  drawbacks  of  the  discrete  (time)  wavelet 
transform  is  the  shift  variance,  since  the  wavelet  coefficients 
of  a  signal  and  a  shifted  replica  of  itself  can  be  very 
different.  A  Euclidean  distance  measure  (introduced  by  Aware, 
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Inc.,  1992)  is  used  as  part  of  the  technique  to  classify  the 
signals  at  scales  where  the  effects  of  shift  variance  can  be 
tolerated. 

Mallat  and  Zhong  (1989)  demonstrated  that  the  maxima 
extracted  from  the  modulus  of  the  wavelet  coefficients  can  be 
used  to  reconstruct  the  input  signal.  That  is,  the  maxima  of  the 
modulus  of  the  wavelet  coefficients  contain  approximately  the 
same  amount  of  information  as  the  original  signal.  Consequently, 
signal  analysis  can  be  performed  based  on  the  wavelet  extrema, 
which  form  a  reduced  signal  representation. 

Thus,  wavelet  coefficients  at  each  scale  are  replaced  by 
their  extrema.  The  reduced  set  is  only  non  zero  where  the  scale 
has  an  extrema,  and  is  equal  to  the  original  value  at  these 
locations.  Wavelet  scale  coefficients,  which  are  not  extrema,  are 
set  to  zero. 


2.  Ranking/Pairing  Algorithm 

The  first  step,  in  computing  a  distance  measure  based  on 
pairing,  is  to  rank  the  peaks  of  the  sets  to  be  compared.  That 
is,  the  peaks  at  each  scale  are  ranked  by  their  amplitudes  (i.e. 
same  procedure  as  is  used  in  median  filtering)  .  Therefore  the 
ranked  sequences  are  ordered  (ranked)  starting  with  the  smallest 
value  and  ending  with  the  largest  value. 

The  second  step  is  to  pair  the  ranked  peaks  of  the  two 
sequences  to  be  compared.  The  peak  with  the  highest  rank  (largest 
amplitude)  in  one  set  is  paired  to  the  peak  with  the  highest  rank 
(largest  amplitude)  in  the  other  set.  The  next  in  rank  (order, 
sequence)  is  paired  to  that  next  in  rank  of  the  other  set.  The 
pairing  scheme  does  not  require  the  number  of  elements  in  both 
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sets  to  be  the  same.  When  a  peak  in  one  set  does  not  have  a 
corresponding  peak  in  the  other,  a  zero  is  inserted  into  the  set 
that  has  a  smaller  number  of  peaks,  and  the  pairs  are  formed  by 
matching  the  remaining  peaks  and  the  zeros. 


3 .  Distance  Measure 

The  third  step  is  to  corrpute  a  distance  measure  for  the 
matched  pairs.  Several  distance  measures  were  initially  tried. 
The  best  results  were  obtained  by  using  the  ranking  and  pairing 
technique.  It  is  noted  that  the  ranking  approach  will  be 
sensitive  to  the  additive  noise  component.  The  distance  assigned 
to  the  pair  is  the  sum  over  the  Euclidean  distances  in  each 
scale.  Thus,  we  compute 


d(al  ,bi) 


(a**  - 


m' 


(3.2) 


where  (k,m)  denotes  the  locations  of  matched  peaks,  aj  and  IP  are 
the  wavelet  extrema  at  scale  j,  and  ajk  and  tPm  are  the  values  aj, 
h  at  temporal  locations  k  and  m,  respectively.  W\/m  is  the 
weighting  factor  at  scale  j  for  the  relative  distances  between 
the  corresponding  coordinates  of  matched  peaks.  The  weighting 
factor,  W\m  ,  is  defined  as 


;  if  k  *  m 
;  if  k  =  m 


(3.3) 


where  njk  and  njm  are  the  coordinates  of  ajk  and  bPm,  respectively. 
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The  similarity  between  two  sets  is  described  in  terms  of  the 
sum  of  the  Euclidean  distances  of  amplitudes  weighted  by  the 
square  root  of  the  relative  distances  between  the  corresponding 
coordinates.  Basically,  a  small  value  |njk  -  n3J  and  |  ajk  -  bPj  , 
implies  a  high  degree  of  similarity  between  a3  and  b3 . 

The  penalty  weight,  W3  k  m  ,  is  determined  by  the  separation 
of  the  peaks  and  has  a  lower  bound  of  unity.  Large  separation 
between  the  matched  peaks  corresponds  to  large  penalty  factors. 
The  distance  measure  also  depends  on  the  amplitude  difference  of 
the  matched  peaks.  The  distance  measure  is  directly  proportional 
to  the  difference  in  amplitude.  It  should  be  noted  that  even 
identical  scale  outputs  with  an  offset  in  time  (due  to  signal 
delay)  would  have  a  non  zero  weighting  factor.  This  problem  is 
eliminated  by  removing  the  offset  time  prior  to  forming  the 
distance  measure.  One  can  do  this  by  correlating  the  signals  to 
be  identified,  with  the  template  and  subtracting  off  the  amount 
of  misalignment  (i.e.  compensate  for  the  amount  of  lag  indicated 
by  the  cross  correlation)  .  A  second  way,  which  was  used  in  the 
work  presented  here,  is  to  line  up  the  dominant  peaks  of  the 
signal  of  interest  and  the  template  in  the  scale  under 
consideration.  This  naturally  lines  up  the  coefficients  (i.e. 
extrema)  in  a  way  that  compensates  for  relative  time  off  set 
between  the  two  sequences  under  observation. 


4 .  Implementation 

There  is  a  total  of  four  different  signal  sets.  Each  set  is 
generated  by  a  different  transmitter.  The  method  outlined  so  far 
was  employed  to  determine  the  transmitting  source  of  the  signal. 
For  this  purpose,  a  tenplate  for  each  transmitter  set  is  needed 
to  measure  the  degree  of  similarity  with  any  signal  of  interest. 
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The  distance  measure  algorithm  is  applied  to  a  tenplate  and  a 
signal  of  interest  at  the  appropriate  wavelet  scale  after 
conpensation  of  the  time  offset. 

The  first  signal  from  each  of  the  four  sets  was  chosen  as  a 
template.  A  small  distance  measure  iirplies  that  the  signal  under 
consideration  comes  from  the  set  that  the  tenplate  represents. 

The  fourth,  eighth,  and  sixteenth  order  Daubechies  wavelet 
functions  were  used  to  conpute  the  WT  of  the  pre-processed 
signals.  The  eighth  order  Daubechies  (Dau-8)  wavelet  functions 
gave  satisfactory  results,  whereas  Dau-4  and  Dau-16  did  not.  The 
Dau-8  wavelet  function,  which  was  used  as  the  mother  wavelet,  is 
shown  in  Figure  5. 


Figure  5.  Eighth  Order  Daubechies  Wavelet  Function 

For  illustration  purposes,  the  distance  values  at  Scale  11 
between  the  first  four  signals  from  Transmitter  1  and  the 
tenplates  are  tabulated  in  Table  4.1.  As  expected,  the  distance 
value  with  Tenplate  1  is  significantly  less  than  that  with  the 
other  tenplates,  since  Tenplate  1  is  the  tenplate  for  the  signals 
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from  Transmitter  1,  and  small  distance  values  imply  high 
similarity. 


signal  1 
signal  2 
signal  3 
signal  4 


Tempi .  1 

0.00 

3.05 

2.61 

1.08 


Tempi .  2 

44.30 

47.23 

40.58 

42.90 


Tempi .  3 

14.77 

16.34 

13.79 

15.48 


Tempi .  4 

107.84 

103.85 
99.78 

105.30 


Table  4.1  Distance  Measures  'd'  at  Scale  11  for  Signals  1-4  of 

Transmitter  1 


As  can  be  seen,  the  distance  values  '  d '  between  signals  from 
Transmitter  1  and  Template  1  are  well  separated  from  those  of  the 
other  templates.  Zero  value  in  distance  indicates  a  perfect  match 
which  occurs  between  Signal  1  and  Template  1.  This  is  expected 
since  Signal  1  was  chosen  as  a  template  for  the  set.  Similar 
results  are  obtained  when  comparing  the  different  signal  sets 
with  the  four  templates  [Payal,1995] . 

A  thresholding  technique  can  be  introduced  to  automate  the 
identification  procedure  by  defining  a  threshold  level  for  each 
tenplate.  For  example,  one  can  compare  the  maximum  value  of  Set  1 
with  the  minimum  values  of  other  templates.  Determining  the 
threshold  levels  is  not  addressed  in  this  study;  however, 
comparing  the  minimum  and  maximum  values  of  the  distance  measure 
is  quite  useful  in  evaluating  the  performance  of  the 
identification  scheme.  The  maximum  and  minimum  values  quoted  in 
the  remainder  of  this  chapter  are  obtained  by  using  all  template 
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and  all  data  sets.  At  Scale  11  with  Template  1,  the  maximum 
distance  value  is  5.3971,  whereas  the  nearest  minimum  distance 
value  of  the  other  three  templates  is  12.52.  This  minimum  occurs 
with  Template  3,  which  indicates  that  Transmitter  1  and 
Transmitter  3  have  a  somewhat  similar  transient  response  (i.e., 
turn-on  behavior)  .  The  ratio  of  the  minimum  Template  3  output  to 
the  maximum  Template  1  output  is  2.319,  which  shows  how  well  the 
Template  1  results  are  separated  from  the  other  template  results. 
This  separation  ratio  is  2.381  and  1.66  for  Scale  10  and  Scale  9, 
respectively.  Lower  scales  (lower  than  9)  were  not  useful  in 
terms  of  the  similarity  measurement,  and  are  not  included. 

Similarly,  the  application  of  the  distance  measure  to  the 
other  three  signal  sets  resulted  in  the  following  separation 
ratios.  Again  these  separation  values  are  the  ratio  of  the 
smallest  wrong  set  distance  and  the  largest  correct  set  distance. 
They  show  how  well  the  signals  of  a  particular  set  are  separated 
from  the  other  sets  under  worst  case  conditions.  The  separation 
ratios  are  3.56,  2.76,  and  1.88  for  signals  from  Transmitter  2  at 
Scales  11,  10,  and  9,  respectively.  The  following  separation 
ratios  were  obtained  at  Scales  11,  10,  and  9:  4.64,  1.98,  and 
3.73  for  signals  from  Transmitter  3;  and  1.75,  1.4,  1.24  for 
signals  from  Transmitter  4.  Signal  versus  template  (distance 
measure)  output  plots  for  all  signal  sets  are  given  in  Appendix 
A.  One  can  obtain  minima,  maxima  ,  and  the  spread  of  the  distance 
measures  very  easily  from  the  plots.  It  should  be  noted  that,  for 
Transmitters  1,  2,  and  3,  nine  signals  were  used  while,  for 
Transmitter  4,  only  six  signals  were  usable. 


5.  Probability  of  Identification 

Signal -to-Noise  Ratios  (SNR's)  for  all  signals  were 
estimated.  It  is  clear  from  Figure  1  that  the  signal  waveform  can 
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be  partitioned  into  three  regions .  The  off  region  is  where  there 
is  no  signal;  the  transition  region  is  where  there  is  a  build-up 
from  off  to  on  state;  and  the  on  region  is  where  the  signal  is  at 
steady  state.  We  can  assume  that  the  off  region  consists  of  noise 
on-ly/  allowing  the  computation  of  the  noise  power.  Noise  and 
signal  coexist  in  the  on  region.  The  signal  power  can  be  computed 
by  finding  the  power  in  the  on  region  and  subtracting  the  noise 
power.  Hence,  we  can  compute  the  SNR's  for  all  sets.  An  average 
SNR  value  of  32.38  dB,  40.87  dB,  38.58  dB,  and  31  dB  was  computed 
for  Set  1,  Set  2,  Set  3,  and  Set  4,  respectively. 

Gaussian  white  noise  is  added  to  all  signal  sets  to  decrease 
the  SNR  levels.  To  obtain  some  statistical  significance  the 
experiments  are  repeated  eleven  times  for  each  SNR  value,  and  the 
probability  of  identification  (Pz)  versus  SNR  is  computed.  The 
Pi,  the  probability  of  identification  is  defined  as  the  ratio  of 
the  number  of  correct  identification  and  the  total  number  of 
experiments . 

The  results  are  shown  in  Figure  6.  If  the  signal  of  interest 
is  from  Transmitter  1  and  its  SNR  is  better  than  16  dB,  then  its 
identification  probability  is  high.  The  minimum  SNR  values 
required  for  a  reliable  classification  at  Scale  11  are  30,  23, 
and  30  dB  for  signals  from  Transmitter  2,  3,  and  4,  respectively. 
Several  line  up  (pairing)  procedures  were  tried  on  an 
experimental  basis,  to  see  if  the  identification  procedure  can  be 
made  more  robust.  So  far  only  the  distance  measure  that  uses  all 
extrema  and  pairs  the  ranked  peaks  seems  to  be  effective  in 
classifying  the  four  signal  sets. 
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Figure  6.  Probability  of  identification  (Px)  of  the  signals 
(a)  Pi  of  signals  from  Transmitter  1  (b)  Px  of  signals  from 
Transmitter  2  (c)  Px  of  signals  from  Transmitter  3  (d)  Px  of 
signals  from  Transmitter  4. 


V.  CONCLUSION 

1.  CLASSIFICATION  OF  PUSH-TO-TALK  COMMUNICATION 

The  main  objective  of  this  study  is  to  determine  a  robust 
wavelet  based  algorithm  designed  to  extract  features  to  identify 
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push- to -talk  transmitters.  Robustness,  compact  signal 
representation  capability,  and  low  computational  complexity  are 
the  main  advantages  of  the  wavelet  analysis,  which  is  used  in  the 
feature  extraction  (identification)  . 

Push-to-talk  communication  recordings  provided  for  this 
research  have  a  common  feature:  They  all  include  a  transition 
from  the  off-to-on  state,  as  well  as  the  on-to-off  state.  The 
turn-on  transition  phase  is  unique  for  each  transmitter,  and  can 
effectively  be  used  for  classification  purposes.  The  recordings 
differ  from  each  other  by  the  waveform  in  the  length  of  the 
transition  region.  The  feature  for  the  classification  of  the 
signals  is  contained  in  the  signal  envelopes .  The  turn-on  part  of 
the  envelope  is  a  transient  and,  hence,  a  wideband  signal. 

The  original  push-to-talk  transmitter  recordings  are  not  in 
appropriate  form  for  WT  analysis.  Thus,  the  recordings  are  pre- 
processed  to  be  usable  for  WT  processing.  Differentiation  of  the 
envelope  of  the  signals  is  used  to  transform  the  data  into 
pulse- shaped  transients.  Filtering  is  applied  to  both  the 
envelope  and  the  differential  to  improve  the  signal-to-noise 
ratio . 

A  distance  algorithm  is  introduced  in  this  work.  It  is  based 
on  an  Euclidean  distance  measure  between  the  wavelet  coefficients 
of  two  data  set  in  terms  of  magnitude  and  relative  position  on  a 
given  scale.  Decisions  about  the  origin  of  the  signal  are  made 
according  to  the  distance  measures  between  the  signals  and  the 
templates,  where  each  template  represents  a  different 
transmitter.  A  small  distance  value  implies  that  the  signal 
belongs  to  the  same  set  as  that  particular  template.  In  its 
current  form,  the  classification  assumes  that  any  signal  of 
interest  belongs  to  one  of  the  four  sets. 

The  distance  algorithm  was  applied  to  four  different  signal 


19 


sets.  The  first  recording  in  any  of  the  sets  is  designated  to  be 
the  template.  Instead  of  using  all  the  wavelet  coefficients,  just 
the  local  extrema  are  used.  Using  only  the  local  extrema  reduces 
the  corrputational  complexity  of  the  algorithm. 

The  distance  measure  between  the  signal  and  the  templates  is 
the  sum  of  all  Euclidean  distances  between  paired  peaks  of  the 
signal  and  the  templates.  It  also  includes  a  penalty  factor  due 
to  the  relative  square  root  distance  between  the  matched  ranked 
pairs  and  their  difference  in  amplitude.  Matching  signal  peaks  to 
the  template  is  performed  by  ranking.  Before  pairing  the  peaks, 
the  maximum  peak  of  the  signal  is  aligned  with  the  maximum  of  the 
templates. This  tends  to  reduce  the  penalty  weight  for  like 
signals  which  are  not  aligned  in  time. 

The  distance  algorithm  described  in  chapter  4,  the  one  that 
uses  all  the  extrema  and  the  pairing  of  ranked  peaks  allows 
robust  identification  of  the  signal  sets.  A  printout  for  all 
programs  is  provided  in  Appendix  B. 


2.  RECOMMENDATION  FOR  FUTURE  STUDY 

The  distance  algorithm  introduced  in  Chapter  4  gave 
promising  results  in  classifying  the  four  signal  sets  provided 
for  this  study.  Three  important  issues  have  not  been  addressed  in 
this  research:  i)  the  template  selection,  ii)  a  threshold 
technique,  and  iii)  incorporation  of  information  from  other 
scales.  In  this  work,  templates  are  chosen  arbitrarily  from  the 
signal  sets.  When  the  signals  to  be  identified  are  from  the  four 
sets,  the  algorithm  is  capable  of  classifying  the  signals. 
However,  if  the  signals  do  not  belong  to  these  sets,  the 
algorithm  will  compute  a  distance  to  each  of  the  templates,  which 
could  lead  to  misinterpretation  hence  misclassification. 
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A  thresholding  technique  and  robust  template  selection 
should  be  the  subject  of  further  study. 

Also,  no  attempt  was  made  to  combine  information  from 
several  scales.  Identification  was  obtained  by  using  just  one 
scale  (i.e.  highest  frequency  location)  .  Typically,  Scale  11 
worked  best.  If  the  distance  information  from  Scales  9  and  10 
could  be  used,  a  potentially  more  robust  identification 
performance  should  be  realized. 
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Appendix  A 

Appendix  A  consist  of  2  sets  of  4  figures  (Figure  A.l  - 
A. 8) .  The  first  four  figures  use  data  at  their  original  SNR  value 
while  the  last  four  figures  use  data  at  an  SNR  level  10  dB  below 
the  original  values.  The  absolute  SNR  levels  of  the  originals  are 
32.38,  40.87,  38.58,  and  31  dB  for  transmitters  1,2,3  and  4, 
respectively.  Figure  A.l  and  A. 5  use  signals  from  transmitter  1. 
Figure  A.  2  and  A.  6  use  signals  from  transmitter  2.  Figure  A.  3  and 
A.  7  use  signals  from  transmitter  3.  Figure  A. 4  and  A.  8  use 
signals  from  transmitter  4. 

These  figures  allow  easy  determination  of  the  maximum  and 
minimum  values  as  well  as  the  mean  and  variability  behavior  of 
each  test  SNR  level.  One  can  fairly  easily  establish  the 
sensitivity  of  each  set  relative  to  its  template  and  to  members 
of  the  other  sets.  The  second  set  (Fig.A.5-A.8)  provide  a  sense 
of  the  degradation  of  the  identification  procedure  as  the  SNR  is 
lowered  by  10  dB.  More  information  in  terms  of  identification 
ability  and  SNR  levels  are  provided  in  the  main  body  (see  Fig. 

6)  . 
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Figure  A.  1  Distance  measures  at  the  output  of  templates  when,  signals  from  Trans¬ 
mitter  1  are  the  inputs.  Horizontal  axis  shows  the  number  of  the  signals.  Each  column 
represents  a  template.  First  row  is  for  Scale  11;  second  row  is  for  Scale  10;  third  row 
is  for  Scale  9. 
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Figure  A.  2  Distance  measures  at  the  output  of  templates  when  signals  from  Trans¬ 
mitter  2  are  the  inputs.  Horizontal  axis  shows  the  number  of  the  signals.  Each  column 
represents  a  template.  First  row  is  for  Scale  11;  second  row  is  for  Scale  10;  third  row 
is  for  Scale  9. 
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Figure  A.  3  Distance  measures  at  the  output  of  templates  when  signals  from  Trans¬ 
mitter  3  are  the  inputs.  Horizontal  axis  shows  the  number  of  the  signals.  Bach  column 
represents  a  template.  First  row  is  for  Scale  11;  second  row  is  for  Scale  10;  third  row 
is  for  Scale  9. 
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Figure  A.  4  Distance  measures  when 
Horizontal  axis  shows  the  number  of  the 
First  row  is  for  Scale  11;  second  row  is  f 
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Figure  A.  5  Distance  measures  when  signals  from  Transmitter  1  are  inputs  and 
their  SNR  values  are  10  dB  lower  than  the  original  data  set.  Horizontal  axis  shows 
the  number  of  the  signals.  Each  column  represents  a  template.  First  row  is  for  Scale 
11;  second  row  is  for  Scale  10;  third  row  is  for  Scale  9. 
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Figure  A.  6  Distance  measures  when  signals  from  Transmitter  2  are  inputs  and 
their  SNR  values  are  10  dB  lower  than  the  original  data  set.  Horizontal  axis  shows 
the  number  of  the  signals.  Bach  column  represents  a  template.  First  row  is  for  Scale 
11;  second  row  is  for  Scale  10;  third  row  is  for  Scale  9. 
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Figure  A.  8  Distance  measures  when  signals  from  Transmitter  4  are  inputs  and 
their  SNR  values  are  10  dB  lower  than  the  original  data  set.  Horizontal  axis  shows 
the  number  of  the  signals.  Each  column  represents  a  template.  First  row  is  for  Scale 
11;  second  row  is  for  Scale  10;  third  row  is  for  Scale  9. 
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Appendix  B 


This  appendix  provides  a  listing  of  the  Matlab  code  used  to 
obtain  the  identification  probabilities. 
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WTSIGNAL.M 


%  WTSIGNAL.M 

%  This  program  removes  the  DC,  takes  the  envelope,  and  filters 
%  the  record  of  a  transmitter.  Then  it  takes  the 
%  differential  and  filters  again  to  obtain  a  signal  applicable  to 
%  wavelet  processing. 

x= input ( '  enter  the  signal  name  • );%all  programs  assume  a  data 

%  vector  of  length  N,  which  is  a  power  of  2 

x=x-mean (x) ;  %  removes  DC  component 

y=envelope  (x)  ,-  %  takes  envelope 

my=asmooth (y, 100) ;%  boxcar  averaging  size  100 

dmy=diff (my) ;%  differential  of  the  envelope 

mdmy=asmooth(dmy, 50) ;%boxcar  averaging  size  50 

%  add  on  one  data  point  (zero)  for  the  lost  one  in  dif ferentiatio 
n 

i=length (mdmy) ; 
mdmy (i+1) =0 ; 
clear  x  y  my  dmy  i 
end 


33 


ASMOOTH.M 


ASMOOTH.M 


function  [y]  =  asmooth(x,L) 

%  boxcar  averager  of  length  L,  creates  as  many  data  points 
%  (length  of  y)  as  the  length  of  the  input  vector  x. 

y  =  []  ; 

if  nargin  ~=  2, 

error ( 'avsmooth:  invalid  number  of  input  arguments. . . ' ) ; 

end 

if  min(size(x))  ~=  1, 

error (' avsmooth :  input  argument  must  be  a  lxN  orNxi  vector'); 
end 

x=x ( : ) ; 

ns  =  length(x); 
y=zeros (ns, 1) ; 

%  average 

%For  the  first  L/2  points 
for  k=l:L/2 

y (k, 1) =mean(x(l:L/2+k, 1) ) ; 

end 

%for  in  the  data 
for  k=L/2+l :ns-L/2 

y (k, 1) =mean (x (k-L/2 :k+L/2-l,  1)  )  ; 

end 

%for  going  out  of  the  data 
for  k=ns-L/2+l :ns 

y (k, 1) =mean (x (k-L/2 :ns, 1) ) ; 

end 
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%  ENVELOPE. M 

function  [y;m]  =  envelope (x) 

%  computes  the  envelope  by  taking  the  absolute  value  of 
%  the  Hilbert  transform 

7  =  []  ; 

if  nargin  ~=  1, 

error (' envelope:  only  one  argument  allowed'); 

end 

if  min(size(x))  ~=  1, 

error (' envelope :  input  argument  must  be  a  lxN  orNxl  vector'); 
end 

x=x  ( : )  ; 

y=abs (hilbert (x) ) ; 
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%  MAP .M 

function  A  =  map(f,N) 

%  modified  version  of  mapdn.m.  It  computes  the  amplitudes  at  the 
%  scale  outputs . 

M  =  length (f); 
n  =  round (log (M) /log (2) ) ; 
a  =  wavedn(f,N) ; 
b  (1)  =  a  (1)  ; b  ( 2 )  =  a  (2)  ; 
for  j  =  l:n-l 
for  k  =  l:2Aj 

index  =  2^j+k+N/2-l; 

while  index  >  2A  (j+1)  ,  index  =  index-2xj  ,-end 
b(index)  =  a(2Aj+k); 
end 
end 
a  =  b; 

for  j  =  i:2x (n-1) 

A (1, j )  =a (1)  ; 
end 

for  j  =  2:n+l 

for  k  =  1 : 2* ( j -2) 
for  m  =  1:2* (n-j+1) 

A(j  ,  (k-1) *2a (n-j+l) +m)  =  a (2* ( j -2 ) +k) ; 
end 
end 
end 
A=A; 
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%  WAVEDN.M 

function  a  =  wavedn(f,N) 

% 

o 

o 

M  =  length (f); 
n  =  round (log (M) /log (2) ) ; 
c  =  dcoeffs(N); 
clr  =  fliplr(c) ; 

for  j  =  1:2: N- 1  ,  clr(j)  =  -clr(j)  ;  end 
a  =  f; 

for  k  =  n: -1: 1 
m  =  2^ (k-1) ; 

x  =  to]  ;  y  =  [0]  ; 

for  i  =  l:m 
for  j  =  1 :  N 

k(j  )  =  2*i-2  + j  ; 

while  k( j )  >  2*m  ,  k ( j )  =  k(j)-2*m 
end 

z  =  a  (k)  ; 

[mr,nc]  =  size(z); 
if  nc  >  1  ,  z  =  z '  ;  end 
x(i)  =  c*z ; 
y(i)  =  clr*z; 
end 

x  =  x/2  ;  y  =  y/2  ; 
a ( 1 :  m)  =  x; 
a(m+l:2*m)  =  y; 
end 


;  end 
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%  DCOEFFS.M 

function  c  =  dcoef fs (N) 

o. 

"o 

nm  =  sqrt ( 2 ) ; 
c  =  zeros (1,N) ; 
if  N  ==  2 
C  =  [11]  ; 

end 

if  N  ==  4 

c  =  [ (1+sqrt (3) ) /4  (3+sqrt (3 ) ) /4  (3-sqrt (3) ) /4  (1-sqrt (3) ) /4] ; 
end 

if  N  ==  6 

q  =  sqrt(10);s  =  sqrt(5+2*q); 

c=  [ (1+q+s) /16  (5+q+3*s)/l6  (5-q+s) /8  (5-q-s) /8  (5+q-3*s)/l6  (1 

+q-s) /16] ; 

end 

if  N  ==  8 

c  =  [.3258030428051  ,1.010945715092  . 892200138246 039575026236 , 
-.264507167369, .043616300475, .046503601071,-0.14986989330] ; 
end 

if  N  ==  10 

c  =  [.226418982583,  .853943542705,1.024326944260,  . 195766961347 ,-.  34 
2656715382, - .045601131884, . 109702658642 , - . 008826800109 , - . 017791870 

102. . 004717427938] ; 
end 

if  N  ==  12 

c  =  [.157742432003, .699503814075,1.062263759882, . 445831322930 ,-. 31 
9986598891, - . 183518064060 , .137888092974, .038923209708, - .0446637483 

31. . 000783251152. .006756062363, - . 001523533805] ; 
end 

if  N  ==  14 

c  =  [.110099430746, .560791283626,1.031148491636, . 664372482211, -. 20 
3513822463, - . 316835011281, .100846465010, .114003445160, - .0537824525 
90, - .023439941565, .017749792379, . 000607514996 ,-. 002547904718 , .0005 
00226853] ; 
end 

if  N  ==  16 

c  =  [ .076955622108, .442467247152, .955486150427, .827816532422, - .022 
385735333, - . 401658632 782 , .000668194093, . 182076356847 , - . 02456390104 
6, - .062350206651, .019772159296, . 012368844819 , - . 006887719256 , - . 0005 
54004548,  .  000955229711, -. 000166137261]  ; 
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end 

if  N  ==  18 

C  =  [.053850349589,  .344834303815,  .855349064359,  .929545714366,  .1883 
69549506, -.414751761802, -.136953549025, .210068342279, .043452675461 
, -.095647264120, .000354892813, . 031624165853 ,-. 006679620227 ,-. 00605 
4960574,  .002612967280,  .  000325814672 ,-. 000356329759 ,  .000055645514]  ; 
end 

if  N  ==  20 

C  =  [.037717157593,  .266122182794,  .745575071487,  .973628110734,  .3976 
37741770, -  .353336201794, - .277109878720,  .180127448534,  .131602987102 
,-. 100966571196, -.04 1659248088,. 046969814 097, .005100436968, - .01517 
9002335 , .001973325365 ,.0028 17686590, -.000969947840, -.0001 64 709006, 
.000132354366, - .000018758416] ; 
end 
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%  DISTANCE.  M 

%  set  up  to  work  with  4  templates  (i.e.  4  transmitters) 
sl=input ( 1  enter  the  signal  name');%  wavelet  transform  output 


kkl = input (' enter  WT  matrix  of  template  1 
%mml= input  (' enter  WT  matrix  of  template  2 
%ttl=input ( ' enter  WT  matrix  of  template  3  : ' ) ; 

%wl=input  ( '  enter  WT  matrix  of  template  4  :  ' )  ; 


%  assumes  scales  6,7,8,9,10,11  are  of  interest 

%(i.e.  data  length  =4096  =  2*12  therefore  11  scales 

%with  only  the  last  6  of  interest  to  this  study 

%  scales  6  to  11  ;  n=scale  #;  highest  scale  =  highest  freq  band 

for  n=6 : 11 

t=length(kkl (1, : ) ) / (2^n) ; 

col=n+2;  %  corresponding  row  for  the  scale 
a=kkl (col , 1 : t : length (kkl (1, : ) ) ) ; 

%b=mml ( col , 1 : t : length (mml ( 1 ,  : ) ) )  ; 

%c=ttl (col, 1 : t : length (ttl (l, :  )  )  )  ; 

%f  =wl  (col ,  1 :  t :  length  (wl  (1,  : )  )  )  ; 

s=sl (col, 1 : t : length (si (1, :))); %  signal  to  be  tested 

%  find  local  extrema 
a=localext (a) ; 

%b=localext (b) ; 

%c=localext (c) ; 

%f=localext(f) ; 
d=localext (s) ; 

%sort  extrema  in  ascending  order 
[tempi, i] =sort (a) ;%  template  1 
% [temp2, j] =sort (b) ;%  template  2 
% [temp3,k] =sort (c) ;%  template  3 
% [temp4,g] =sort (f ) ;%  template  4 
[x,m]  =sort  (d)  ; 

%  Difference  measures 

shiftl=m(length(x) ) -i (length(a) ) ;%shift  for  line  up  of  max  peaks 
%shift2=m (length (x) ) -j (length (b) ) ;%shift  for  line  up  of  max  peaks 

%shift3=m (length (x) ) -k(length(c) ) ;%shift  for  line  up  of  max  peaks 

%shift4=m (length (x) ) -q(length(f ) ) ;%shift  for  line  up  of  max  peaks 

wl= [abs (m-i-shif tl) ] ;  %penalty  weights 
%w2=  [abs  (m- j -shift2)  ]  ,•  %penalty  weights 


40 


D I STANCE. M 


%w3= [abs (m-k-shif t3) ] ;  %penalty  weights 
%w4= [abs (m-q-shif t4) ] ;  %penalty  weights 

wl (find(wl==0) ) =ones (1, length (find (wl==0) ) ) ;%modific.  for  no  0 
%w2 (find (w2==0) ) =ones (1 , length (find (w2==0) ) ) ; %modif ic .  for  no  0 
%w3 (find (w3==0) ) =ones (1, length (find (w3==0 ) ) ) ; %modif ic .  for  no  0 

%w4 (find (w4==0) )=ones(l, length (find (w4==0) )) ;%modific.  for  no  0 

%  the  smaller  the  distance  dl,  d2,  etc.  the  more  likely  the  test 
%signal  belongs  to  the  transmiter  1,2  ,  etc. 
dl (n) =sum(sqrt (wl .* ( tempi -x) . A2)  )  ; 

%d2 (n) =sum(sqrt (w2 . * (temp2-x) .^2) ) ; 

%d3 (n) =sum (sqrt (w3 .* (temp3-x) . ^2) )  ; 

%d4 (n) =sum(sqrt (w4 . * (temp4-x) .A2)); 

end 

end 
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%  LOCALEXT. M 

function  [ou, k] =localext (sg) 

%  LOCALEXT  extracts  the  local  extrema  of  a  vector 
%  nonextrema  are  set  to  zero 

%  [Y,K]  =localext (x)  will  return  the  number  of  deleted  samples 

o, 

o - —  —  —  —  — - - — - —  — - — - 

%  copyright  1994  by  Universidad  de  Vigo 

%  under  GNU  conditions 

%  Author:  Sergio  J.  Garcia  Galan 
%  e-mail :  Uvi_Wave@sc .uvigo . es 

o, 

o  - — - 

l=length(sg) ; 
sg=sg ( : ) ' ; 
ou=zeros  (1,1) 
ax=abs (sg) ; 

MX=max  (ax)  ; 

MX=MX+1 ; 

OU ( 1 ) =sg ( 1 ) ; 
ou  ( 1 )  =sg  (1)  ; 
k=0 ; 

for  i=2 : 1-1 

if (sg(i) >sg (i-1) ) &(sg(i) >sg(i+l) ) 
ou (i) =sg (i) ; 
k=k+l ; 
end 

if (sg (i) <sg(i-l) ) &(sg(i) <sg(i+l) ) 
ou (i) =sg  (i)  ; 
k=k+l ; 
end 

end 

k=l-k; 
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